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A new numerical method is introduced for calculating the polarizability of an arbitrary dielectric 
object with position dependent complex permittivity. Three separate numerical approaches are 
provided to calculate the dipole moment of a nanoparticle embedded in a dielectric matrix in the 
presence of an applied electric field. Numerical tests confirm the accuracy of this method when 
applied to several cases for which an exact solution is available. This method is especially well suited 
for the calculation of absolute Raman scattering intensities due to acoustic phonons in metallic and 
dielectric nanoparticles embedded in transparent matrices. 
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I. INTRODUCTION 

Spectral features of low frequency Raman scattering 
from nanoparticles (NP) can be explained in terms of 
the vibrational frequencies of acoustic phonons which are 
confined in those NPi Calculation of the displacement 
fields of the modes permits prediction of the Raman se- 
lection rules. Variation of Raman intensity when parallel 
or perpendicular polarization is selected can also be un- 
derstood. 

Past theoretical work on the intensity of Raman scat- 
tering from NP 2 ' 3 ' 4 has not produced quantitative esti- 
mates of scattering intensities. However, key qualitative 
features, such as the distinction between NP volume and 
surface mechanisms^ have been pointed out. In the NP 
volume mechanism, deformation potential coupling mod- 
ulates the bulk dielectric response^ In the surface cou- 
pling mechanism, changes of the NP's size or shape mod- 
ulate the NP's polarizabilityiS 

Calculation of the absolute intensity of low frequency 
Raman scattering from an NP requires the values of the 
polarizability matrix, a, of the NP, in particular the mod- 
ulation of the polarizability with time due to the acoustic 
phonon degrees of freedom. Our approach incorporates 
both the NP volume and surface coupling mechanisms. 
In addition, our general approach includes an additional 
qualitative mechanism which has not been previously 
pointed out: Variations of the density of the matrix ma- 
terial surrounding the NP will lead, through deformation 
potential coupling, to a new matrix volume mechanism. 

It is a reasonable approximation in many cases to 
hope to determine the polarizability of an object in 
terms of a local permittivity function, e(r) such that 
D(r)=e(r)E(r). Our discussion only applies to this case. 

In some kinds of NPs this approach will not work. It 



is not always the case that the dielectric response of an 
object can be described in terms of a local permittivity 
function. An example of this is when the Raman scatter- 
ing of a NP is dominated by acoustic phonon modulation 
of electron hole excitons, as in semiconductor NPs when 
the Raman laser is close to the exciton energies. In ad- 
dition, local dielectric response cannot be expected for 
distance scales comparable to the Thomas-Fermi length. 
Such situations will not be considered further in what 
follows. 

What is needed is a method which can handle: (1) 
small variations in shape of an object; (2) complex per- 
mittivity; (3) permittivity which is an arbitrary function 
of position; (4) non-spherical shape; (5) permittivity with 
a large magnitude. 

Quite a number of methods are available to calcu- 
late the polarizability of an object. However, none of 
these methods are suitable for the requirements of low 
frequency Raman scattering from NPs. 

The exact solution due to Mie& is for homogeneous 
spherical objects only. The most general method is Fi- 
nite Difference Time Domain (FDTD)ii However, FDTD 
cannot handle very small changes in the object since it 
must use a coarse grid of spatial points. The discrete 
dipole approximation (DDA)2*2iiS requires a mesh to ap- 
proximate the object, and is not suitable to reflect small 
changes in shape due to vibration. The dipole-induced- 
dipole (DID) method 3 - is applicable only when the per- 
mittivity is small. 

No presently available numerical method meets all of 
these criteria. This paper introduces a new method which 
satisfies all of these requirements. The application of this 
method will permit quantitative estimates of the inten- 
sity of low frequency Raman scattering due to acoustic 
phonons in NPs. 

This paper is restricted only to the problem of calcu- 
lating the polarizability tensor for a given static configu- 
ration of a NP. Repetition of this method for a sequence 
of configurations associated with the motion of an acous- 



2 



tic phonon will allow the modulation of the polarizability 
tensor to be determined. This leads directly to the scat- 
tered Raman intensity, which is the ultimate justification 
for this work. 

NPs are very small compared to the wavelength of the 
laser light used to excite them in Raman and Brillouin 
light scattering experiments. Thus, at an instant of time 
the electric field in the region enclosing a NP may be 
regarded as a uniform static electric field. The NP has 
a dielectric constant which differs from that of the sur- 
rounding glass matrix. For this reason there will be a 
spatial variation of the electric field in the interior and 
vicinity of the NP. This electric field induces polariza- 
tion and consequently bound charge on the NP which 
leads to a dipole moment. It is this electric dipole that 
oscillates so as to radiate, emitting Rayleigh scattered 
light (with the same wavelength as the incident light). 
Acoustic vibrations of the NP lead to variations in the 
dielectric response of the NP. Their frequency is very low 
compared to the laser. As the NP slowly vibrates, its 
dielectric response leads to Raman or Brillouin scattered 
light whose frequency is slightly above or below that of 
the incident laser. 

An electrostatic (also called "quasistatic" ) description 
of a system is justifiable if the characteristic frequencies 
applied are much less than the speed of light divided by 
the diameter of the system. For a NP, this intrinsic fre- 
quency would be roughly 10 17 Hz, whereas the frequency 
of light is below 10 15 Hz. Consideration of retardation 
effects is not important if we are considering small NPs. 
In addition, note that the NP oscillates very slowly com- 
pared to the incident electric field. 

To help understand the physical situation in question, 
consider a single instant in time. Consider also a re- 
gion surrounding the NP over which the incident elec- 
tric field is approximately constant. The incident electric 
field E mc has components E incx , E incy , and E mcz . If the 
region under consideration contained only glass without 
a NP, then the electric field would be constant within 
the region. The effect of the variation of the dielectric 
response of the NP material relative to that of the glass 
is that the electric field varies in and close to the NP. 

In the following sections, a numerical method is intro- 
duced for determining the resulting electric field associ- 
ated with the NP as well as the dipole moment of the NP. 
This can be used to find the polarizability of the NP. 



II. DIPOLE RADIATION 

Before continuing on to the main point of this paper, 
which is the calculation of the polarizability of a NP, we 
will here present the key relationships that will permit 
these results to be applied to the problem of the calcula- 
tion of the intensity of Raman scattering. 

First, note that, for a small object, the radiation emit- 
ted by it is completely dominated by its dipole moment, 
and that higher order moments such as quadrupole are 



negligible. If an object in vacuum has dipole moment 
z-component varying with time as p z = p oz cos(wi) then 
the radiated power per unit area of the detector isii 
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where p a and c are the permeability and speed of light, 
respectively, of free space. r pc i is the distance from the 
NP to the photon detector. is the angle between the 
axis of the dipole and the ray from the dipole to the 
photon detector. 

Even though the results of this paper may be applied 
to the situation of a NP in vacuum, it is more common 
in experimental situations that the NP is embedded in a 
macroscopic matrix such as a block of glass. It is straight- 
forward to adapt Eq. (JTJ to this situation: 
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where p m is the permeability of the matrix and c rn is 
the speed of light in the matrix. c m = 1/ '-^/jMn^m an d 
e m = e mro e 0l where e mro is the relative permittivity of 
the matrix. e is the permittivity of free space. 

As an aside, we add the "o" to e mro to signify that this 
is the equilibrium value of this quantity, which is neces- 
sary because our later calculations can involve situations 
where the vibrations of the NP cause variations in the 
density of the matrix material and consequently result in 
changes of the permittivity of the matrix in the imme- 
diate vicinity of the NP. The formalism in later sections 
allows the permittivity of the matrix in the immediate 
vicinity of the NP to be a function of position. 

When dealing with the macroscopic behavior of di- 
electrics, charge may be viewed as "free" , "bound" , or 
"total". Charge which is artificially added to a preex- 
isting neutral dielectric material is certainly "free" . The 
response of the dielectric is to rearrange its own charge 
so as to partly screen the "free" charge. Concentrations 
of this redistributed charge are "bound" charge. Specif- 
ically, if pf, pb, and p are, respectively, the free, bound, 
and total charge densities, then pf +pb = p, V- E = p/e a 
and V • D = pf. 

A point charge qj ( "f ' stands for "free" ) in vacuum cre- 
ates an electric field with magnitude E = kqf/r 2 where 
k = l/(47re ). If the same point charge qf is artificially 
added into an initially neutral block of dielectric material 
of relative permittivity e mro , it creates an electric field 
E = kqf /(e mro r 2 ). The reduction in E is due to bound 
charge qb = —{{f-mro — 1) / '^mro)qf ■ The total charge con- 
tained in the vicinity immediately surrounding the free 
charge is q = q f + q b , and is given by q = qf/e mro . 

In like manner, when speaking of a point dipole in a 
dielectric matrix, we must distinguish between the free, 
bound, and total dipole moments. Let these three (vec- 
tor) quantities be denoted respectively by p/, pf,, and 
p. For specificity, suppose that the dipole is oriented 
along the z-axis with respective moments pf z , pb z , and 
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p z . If a free dipole p f z is artificially inserted into a neutral 
dielectric, it will be screened by bound dipole moment 

Pbz — ((^mro 1) /^i7iro)Pfz- 

In subsequent sections of this paper, we will always 
be dealing with situations where there is no free charge 
present. As a result, all of the dipole moments which we 
will calculate involve bound charge only. When these 
dipoles oscillate, they radiate electromagnetic energy. 
This is the fundamental mechanism for all Rayleigh and 
Raman scattering from NPs. It is important to note that 
Eq. and Eq. cannot be used to find the energy ra- 
diated from such NPs because the relative permittivity 
of the dielectric surrounding the NP must be considered. 

The correct way to obtain the radiated energy from a 
NP embedded in a dielectric matrix is as follows. Given 
the bound dipole moment pb z = p z , we find the equiva- 
lent free dipole (p/ 2 ) which would create the same fields 
in the immediate vicinity of the NP. These are related by 

Pfz ^-vnroPz- 

Thus, in the absence of any free dipole moment, for 
an oscillating point dipole p z = p oz cos(ujt) in a dielectric 
matrix, the radiated power density is: 
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III. SPHERICAL HARMONIC TRANSFORM 

Consider a NP that is approximately spherical in shape 
and approximately centered at the origin of the coordi- 
nate system. The electric field is E(r). Because of the 
static approximation, E = — W where V(r) is the elec- 
tric potential. 

The motivation for introducing the word "approxi- 
mately" twice above is to justify the use of a spherical 
harmonic expansion for the potential V which is dom- 
inated by components with slow angular variation. In 
such a situation, we can hope for a useful approximation 
with a finite number of spherical harmonics. However, 
with generality any potential whatsoever could be ex- 
pressed in the form: 



v(r,e, 



at m {r)Si m (6,4>) 



(4) 



where Si m {9, <fi) are "real spherical harmonics" which are 
real- valued functions here defined as: 



Sim — 

Seo = 
S lm = -V2lm(Yl 



[to > 0] 
[to = 0] 
[to < 0] 



The Y["-(9,(j)) are conventionally defined complex val- 
ued spherical harmonic functions^ 
There are orthonormality conditions 



StrnShM Sin 



iLOmM 



(5) 



where &y is the Kronecker delta. In particular, Sqo — 
0.2821, Sio ~ 0.4886 cos (9, S lt ~ 0.4886 sin 9 cos 0, and 
S\ -l ~ 0.4886 sin 9 sin 4>. Finally, note that: 



V 2 S, 



£(£+1) 



Sr. 



(6) 



In macroscopic electrostatics the permittivity e and the 
fields E, D and the polarization P are related through 
D = eE and D = e D E + P so that P = e (e r - 1)E where 
e a is the permitivity of free space and e = e r e , where e r 
is the "relative permittivity." 

The permittivity is a function of position, varying for 
two reasons: (1) The permittivity of the NP is different 
from the permittivity of the surrounding glass matrix. 
(2) Vibrations of the NP will lead to elastic strains which 
will cause small time dependent variations of the permit- 
tivity. (These periodic variations have a time scale on 
the order of 3 ps, which is 1000 times longer than the 
period of the incident laser light beam) 

Motivated by the roughly spherical shape of the NP 
about the origin, a real spherical harmonic expansion is 
employed: 



e(r, 9, (f>) = e }^ qe m (r)Se 



(7) 



It also turns out to be convenient to introduce the 
logarithm (base e) of the relative permittivity b = 
log(e/e )=log(e r ) which has the expansion: 



log (e r (r, 6,4>)) = b(r,i 



(8) 
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The permittivity can, in general be complex- valued, 
so the complex analytic continuation of the logarithm 
function will be used. Note that the coefficients Q m (r) 
and qem(f) are dimensionless, while the ai rn (r) are in 
volts. 

The permittivity is assumed to be initially specified at 
all points within the NP and the matrix. The coefficients 
are found using: 



qLM(r) = — I e(r,9,<p)S LM sin 
e J 



(9) 



and 



c LM (r)= [\og(e r (r,9,<t>))S LM sm6d9d(t> (10) 

While V(r) may have non- analytic variation in the 
vicinity of the surface of the NP due to rapid spatial 
changes in e(r), we suppose that V^(r) is smooth at the 
origin with the following series expansion: 



lim V(r, 



de m r Si m ( 



(11) 



Finally, we suppose that sufficiently far away from the 
NP the permittivity again becomes constant. In this re- 
gion the electric field is supposed to reach a spatially con- 
stant value, but formally the potential can be expanded 
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in the large r limit as: 
lim V(r, 0, q 



(12) 



The x, y, and z incident (far away) electric field 
components are related to eio, en, and ei_i as fol- 
lows: Eincz — -0.4886 eio, E mcx ~ -0.4886 en, and 
E lncv ~ -0.4886 ei_i. 



IV. INTEGRATING THE COULOMB 
EQUATION 

The basis for this paper's calculations is that there are 
no free charges associated with the NP, although there 
will be some bound charge as a result of the induced 
polarization. Thus pf — 0, in which case V • D = but 
D = eE, so V ■ (eE) = 0. Note further that E = - W, 
so that V • (e(W)) = 0. Differentiating, 



eVV + (W) • (Ve) = 



(13) 



Next, divide this through by e(r). If b = log(e/e ) = 
log(e r ) then V6 = (Ve)/e. In this case 



W 2 V + (W) • (Vb) = 



(14) 



The boundary value problem that we are solving is as 
follows: Supposing b(r) to be known, it will be Eq. i|14|) 
which will be solved to yield V(r). The boundary condi- 
tion is that E approaches a specified constant value Ei„ c 
in the limit of large r. 

To do this, we next insert the real spherical harmonic 
expansions of V and b into Eq. (|14f) . The result is mul- 
tiplied by a given Se m (0,<p), and then integrated over 9 
and <fi. This yields a set (indexed by £ and m) of coupled 
ordinary differential equations as follows: 

- r 2 a" m = 2ra' im - £(£ + l)a im 



LAI Xfj, ^ 

K(£m\LM;\fi)\ (15) 



where 



and 



H(£m; LM; Xfx) = / S em S LM S Xfi dn (16) 
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Slm^t-Sx^cM (17) 



<d(j) d<j> 



where dfl denotes sm8d9d(j). The constants H(;;) and 
if (|; ) have to be numerically integrated in advance and 



stored in a lookup table. a" m is the second derivative and 
o! im is the first derivative of ai m (r). 

For the initial small r value in each integration of 
Eq. (|15fl . only one of the constants di m (from Eq. (|11|) ') 
will be 1, and the others are all zero. (Later on, we will 
be able to determine the actual values of all of the di m 
so as to satisfy the large-r boundary conditions of the 
original boundary value problem.) 

For a given choice of initial conditions in terms of the 
constants dp m , these coupled second order ordinary differ- 
ential equations, (Eq. (|15fl ). can be integrated outwards 
from r=0 until large r is reached, at which point the 
constants eg m can be found. 

Because of the linearity of the equations, these are re- 
lated by 



EE* 



ImLMULM 



(18) 



where F is a matrix. All of the coefficients of F can be 
calculated through successive integrations of the coupled 
equations (i.e. varying £ and m). 

Suppose the desired solution corresponds to the case 
where eio — — 2.047 Ei ncz and all other eg m are 0. Let G 
be the inverse of the matrix F . Then 



EEs 



tmLM&LM 



(19) 



The next step is to repeat the integration of the cou- 
pled differential equations with the resulting values of 
dim, in which case all of the values of ai m (r) can be 
found for all values of r. Using Eq. (0J, this then yields 
V(r), from which the electric field, polarization, bound 
charge density and dipole moment can all be calculated. 



V. OBJECTS IN VACUO 

In order to understand how the dipole moment of a 
NP in a dielectric is calculated later in this paper, we 
first review the simpler case of a NP in a vacuum. The 
dipole moment of the NP is p, with Cartesian compo- 
nents p x , p y , and p z . Since there is no free charge, the 
dipole moment arises from the bound charge, whose vol- 
ume density is Pb(r). The polarization of the material 
is P(r), and pb = —V • P. Since polarization is dipole 
moment per unit volume, dipole moment can also be cal- 
culated from: 



P = / Prf 



(20) 



But p can also be obtained directly from its definition, 
given that p = pt in this situation: 



p = J rp b d 3 r 
In addition, since V 2 V — — p/e Q 

p = e / r\7 2 Vd 3 r 



(21) 



(22) 
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Equations i|20|l . (|21|l . and l|22[l will serve as the basis 
of three independent numerical methods for calculating 
the dipole moment (and hence the polarizability tensor) 
of a NP. These methods are presented in Sections IVIII 



IVIIII and IIXI respectively. This threefold redundancy of 
our method serves as a check on the numerical reliability 
of its results. 

To see that Eq. (|2*U|l and Eq. l(2*T|) give the same answer, 
consider an arbitrary closed region U enclosing the NP 
in its interior. Let dU denote the surface of U . In this 
case, P is zero everywhere on dU . 

Next, consider the vector field B = zP, where z is the 
z-component of position. According to the divergence 
theorem: 



f B dA = f (V-B)d 3 r 
Jau Ju 



(23) 



Therefore, /(V-B)<i 3 r = 0, since P = on the boundary 
of this closed region. 



However, since V • (zP) = P z + zV • P thus 



p = / vp b d A r = 
'u 



Prf 3 r 



(24) 



VI. EMBEDDED DIPOLES 

Unfortunately, the formulas of Section [V] do not all 
apply to the situation of interest for an embedded NP, 
where the surrounding material (a glass matrix, for ex- 
ample) has a susceptibility, so that the polarization in 
the matrix is nonzero, and must be taken into account 
when finding the dipole moment of the NP. The dipole 
moment that is relevant to light scattering experiments 
on NPs is the additional dipole moment produced as a 
result of the presence of the NP. The relative permittiv- 
ity of the glass matrix without density variations is e mro . 
The polarization of the glass matrix in the absence of the 
NP would be as follows: Pi nc = e (e mro — l)Ej„ c This is 
uniform, so V • Pi nc = and consequently the formula 
(Eq. (|21|0 p = J rptd 3 r can still be used in this more 
general case. 

Because P is nonzero throughout the matrix, Eq. H2U|) 
cannot be used for a NP embedded in an infinite matrix 
(since the integral would diverge). Here, we derive a 
formula in terms of P that can be used to find the dipole 
moment in this case. 

Once again, let us define the vector field B to be 
B = zP. Then choose a spherical region U of radius R 
centered on the NP, where R is large enough that density 
fluctuations due to vibrations of the NP are negligible, 
and also large enough that fields due to multipole mo- 
ments beyond dipole order can be ignored. Again, we 
can use the divergence theorem, Eq. (|2*5)l . 

First, we evaluate the left side of Eq. - We need to 
know B on the spherical surface dU that surrounds the 
NP. The NP has dipole moment z-component p z . The 



potential due to this dipole (not including that resulting 
from the incident field) is>ii 



V di . 



1 p z COS ( 



(25) 



47re r 2 

The subscript 'dip' stands for 'dipole'. [Note: It may 
be tempting to insert a factor of e rnro into the denomi- 
nator here. But remember that this is not the potential 
due to free charge embedded in a dielectric. There is no 
free charge. To find the potential due to a dipole, we can 
apply the formula V 2 V = —p/e where p is the charge 
density, the sum of free and bound charge.] The electric 
field created in the matrix by the dipole moment of the 
NP has radial component E<n pr = —dVdi P /dr 

1 2p z cos9 

Edipr = -. o (26) 

4ne r J 

The part of the polarization in the matrix created by 
the NP has radial component 



P, 



dipr 



(emro - 1) 2>Pz COS( 

47T r 3 



(27) 



The vector field B is the sum of the part created by 
the NP and the part due to the incident applied field: 



B B^^p ~\~ B^ nc 
The radial component of B^j, is: 

(emro - 1) 2Pz COS 2 I 



4?T 



(28) 



(29) 



B^ip-dA ig equal to Bdi pr dA. We now want to integrate 
this over the sphere of radius R. 



B dip • dA = I B dipr R 2 sv&6d6d<f) 



(emro - 1) 2p z COS 2 6 2 . 



47T 



R 2 



R 2 sin 6d6d(j) 



j \e mro -l) Pz cos 2 OsinOdO 

2 (d-rnro 



-Pz 



As for the incident field part Bi„ c , it is: 



ze (e r . 



(30) 



(31) 



To evaluate the surface integral J Bi nc • dA we need 
only the radial component: 

Bincr = e (e mro - l)r cos 9E incz cos 9 (32) 

where Ei ncz is the z-component of the incident field. We 
need to integrate this over the same sphere of radius R: 

J B inc ■ dA = e (e mro - 1) J E incz R 3 cos 2 9 sin 9d9dcj) 



4tt 



(R e (c 7nr0 l)-^mcz) 



P d 6 r 



(33) 
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where Pi ncz = e (^mro — ^)Ei ncz is the polarization in the 
matrix if the NP were not present. 
We now combine results: 



B dA — J B dip ■ dA + j B mc • dA 

= 2(£mr °~ 1) p z + / P lncz d 3 r 



(34) 

Once again, because V • (zP) — P z + zV • P then 
J B • dA = J (V • B)d 3 r 

= J z(V • P)d 3 r + J P z d 3 r (35) 
Next, since pb = — V • P and p z = J zpi,d 3 r, we get 
2{ - emr l~ l) p, + J P mcz d 3 r = -Pz+J P z d 3 r (36) 



and finally 



Pz 



{Pz - P tncz )d 3 r (37) 



26 mro -\- 1 

which can then be vector generalized as follows: 

P = t ^— f(P - e (e mro - l)E mc )d 3 r (38) 

Z€ mro -j- 1 J 



VII. DIPOLE MOMENT: POLARIZATION 

Finally we address the central problem of calculating 
the dipole moment of an object with inhomogeneous per- 
mittivity. The object has an arbitrary shape. In the case 
of a vibrating NP, the region of inhomogeneous permit- 
tivity extends outside the NP into the glass matrix since 
density variations resulting from the vibrations will affect 
the permittivity of the glass as well as that of the NP. 
Only at sufficient distance from the NP will the relative 
permittivity reach its homogeneous value of e mro . Note 
that e m ro is the square of the index of refraction of the 
glass. 

There is some given incident electric field Ej nc . By 
repeating calculations of p for different orientations of 
Ei nc we can obtain the polarizability tensor a. 

We begin by calculating p z using Eq. IpTFjl. We first 
make the substitution E z = —dV/dz, yielding 



3e D 



-y _ £r )"^7 + € rnro)E incz ^d 3 r 



2e„ v 

(39) 

The following seven steps will re-express Eq. (|3*§1) in 
terms of an explicit quadruple summation over indices £, 
m, L, and M, as shown in Eq. (|44l) : 

Step 1 : Recall the expansion of e r in terms of the sum- 
mation of the qi m in Eq. Q. 



Step 2: We want to put the summation for (1 — e r ) 
into a neat form. Note that 1 = y/inSoo, so that 



oo +1 



(1 - e r ) = [V^k6io - qt m ]S tm (40) 



t=0 m=-l 



Step 3: Note that 



8V n 8V sin9dV 

— = cos 9- — 

az or r 89 



(41) 



Step 4: Recall the expansion of V in terms of the sum- 
mation of the ai m in Eq. ffi. 

Step 5: Take the derivative with respect to r: 



= YJ2 a 'LM(r)S L M(e 

L=a m 

Step 6: And with respect to 9 as well: 



(42) 



L=Q M v 7 

Step 7: Substitute Eq's. gUJ), ©, and |g3J| all into 
Eq. so as to get: 



3e c 



2e r 



ImLM 

qem{r)]Sim(0, <p) cos 9a' LM (r)S LM {9, 4>) 

AttSiq - qg m {r)]Se m (9, (t>)^^-a LM (r) ® S ^„' 

r 89 



(1 - e mro )E mcz y z drdVL (44) 

where d 3 r has been replaced with r 2 drdfl and dfl = 
sm0d0d(j). Next, we want to carry out the angular part 
of the i ntegration. First, note that cos 6* ~ 2.047 Siq = 
y / 47r/3S'io and second note that sin 6* = —(d/d9) cos 9 = 
— An/ 3(8/ 89) S io . Third, refer to the definitions of 
H(; ; ) in Eq. (JTBJ and K(\; ) in Eq. QTfl. Putting these 
all together, 



3e c 



i. 



47T<5£o 

fa LM 

%m ]^ M ff(LM;10;im) 
+ ia iM i^(£m|10;iM) 



+ 47r Seo S m0 5 L1 5mo(1 - e mro )B mC z |r 2 dr (45) 

Note that Eq. 145fl still applies in the situation where 
E incz =0 and E incx or E incy is nonzero. This would be 
the situation when calculating the off-diagonal elements 
of the polarizability, a. 
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The other components p x and p y can be found by anal- 
ogous expressions: x <^> 1 and y — 1. 



3e Q 



2e„ 



4tT<5/ 

I V 3 L 

+ ia L Mi^(^m|ll;LAf)) 
r / 



\r 2 dr (46) 



3e G 



2e r , 



17 fa LM 

- fc ]( a ' iM i/(LM;l-l;H 

+ -a LM K(£m\l-l;LM)) 
r / 



+ 47T ^£0 <5 m <^L1 <>M0 \r z dr (47) 

The dipole moment p and applied electric field E,;„ c 
are related by p = aEi„ c where a is the polarizability 
tensor of the NP. 

It makes things more compact to label vector 
components by the subscript p which runs from 
— 1 to 1. For example, p has components de- 
noted generally by p M and specifically po and 
pi, where the correspondence to usual Cartesian 
notation is as follows: (j>i,p-i,po) = (p x ,p y ,p z ) and 

In this way, it is possible to write a single formula 
which can calculate any of the three components of the 
dipole moment of the NP, where p G {1,-1,0}: 



3e Q 



2e r . 



EE 

lm LM 



47r^£o 



qt m {r)](a' LM {r)H{LM- l^lm) 



+ -a LM (r)K(£m\lp;LM) 



)E lnc Ar 2 dr (48) 



4tt(1 



VIII. DIPOLE MOMENT: CHARGE 

For redundancy, the calculation of Section IV 111 is re- 
peated, but this time based on the charge-based expres- 
sion for dipole moment, Eq. I|21l) . The bound charge 
density is pb = — V • P where the polarization P is given 
by P = e D (e r - 1)E. Also, E = - W. So 



Pb = e V ■ ((e r - 1)W) 



(49) 



Pb = e Q {(e r - 1)VV + (W) • (Ve r )} (50) 
From this, and using Eqs. Q and 



Pb 



= e o j^^fem - V^^mJftmV 2 (a LM S L m) 



■f E E Q^iQtmSem) -Q^{ a LM S ' lm) 
EE -2-in(<UrnSl m )-^(a LM S L M) 



EE 



r z sin 



(qtmSern) -^jia-LM S LAl) } (51) 



Carrying out these derivatives, and also using Eq. |JBJ, 
we get 



Pb 



1 



inSm)^[r 2 a'[ M + 2ra' LM 



- L(L + l)aLM]Sg m SLM 

E E Qtm a LM S e™ S LM 



EE %maiM (^) 



dSg m 8Slm 
~d0 W~ 



+- 



1 dSc m dS. 



LM 



sin 



(52) 



This last equation can now be inserted into Eq. 121|) 
to obtain the dipole moment. 

Suppose we want to calculate a single component of 
the dipole moment, p^ where p could be -1, or 1. If 
p = then we multiply pi, by z — rcos(8) ~ 2.047 vSiq. 
In general we multiply by 2.047 rSi^. The integrand is 
then integrated over all space. It is convenient to carry 
out the angular integrations first: 



EEj (qem-V^5 m )[r 2 a LM +2ra' LM 

im LM 

L(L + l)a LM }H(lp; LM; Im) 
+ r 2 q' lm a' LM H(lp;LM;em) 

+ qe m aLMK(lp\LM;£m)>rdr (53) 



Results obtained using this equation can be compared 
to those obtained from Eq. (|4*5|> . This provides a guard 
against programming errors and numerical problems in 
the integrations over r. 



IX. DIPOLE MOMENT: POTENTIAL 

A third independent method for obtaining dipole 
moment is presented in this section. Starting from 
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Coulomb's law in differential form: V • E = p/e Q , and 
using p = J rpd 3 r as well as E = — VV, we get: 



p = -e J r{V 2 V)r 2 drdn 



(54) 
and real 



Next, use Eq. 10} to express V in terms of < 
spherical harmonics. Also, specialize to p z . 

Pz = -e / ^2 r cos ^ [a-e m Si m ]r 2 drdfl (55) 

em 

Recall that cos ~ 2.047 510. Also, make use of 
Eq. So this becomes: 



Pz = -2.047 e a / y*[— — (T*—at m ) 



_|_ 1) 

~2 — a*m]<SioS«mf' drdfi (56) 

It is now apparent how to generalize from the z-axis to 
the fi axis where \i is -1, 0, or 1. Also, the angular inte- 
gration can be performed (employing Eq. 10), followed 
by the summation over £ and to: 




rV^ + 2rV 1/I 



2rai^]dr 



(57) 



Equation (|57|l for dipole moment is far simpler than 
Eq. (|48|l and Eq. (j53(l . In particular, the qi m coefficients 
are no longer required. In theory the range of integration 
is from r = to oo. However, Eq. Q57f) is more susceptible 
to numerical error buildup in the high r region compared 
to the other two methods, and the upper limit of inte- 
gration has to be carefully restricted. The integration 
step should also be kept small. Even so, it tends to give 
much more accurate results than the other two methods 
methods when tested on prolate ellipsoids. 



the ideal shape. In addition, this smoothing should be 
done in such a way that the derivatives of permittivity 
are also smooth. This is implemented in a CH — h pro- 
gram named varpr27e. However when the eccentricity is 
nonzero the functions become smooth, and in these cases 
even zero-order smoothing (continuous, but 1st deriva- 
tive discontinuous) seems to work. 

The coefficients cg m and qg m that are obtained from 
the smoothed permittivity function are now smooth func- 
tions of r. Even so, they change sharply near certain 
points. It is necessary to use integration with variable 
step size where a much smaller step is used near the dif- 
ficult points. 

It is practical to first obtain the functions cg m and 
qg m . These functions of r can be stored in a look-up 
table. However, because of the sharp variation at certain 
points the step size between entries in the table must be 
variable. This is taken care of in the program varpr27e 
which generates a file as output which is read by varpr28j. 

The lookup table for H(;;) and if(|;) can grow to 
several hundred entries, and access time can be mini- 
mized through use of a hash index (likely first guess) 
method as implemented in the functions H(;;) and K(\; ) 
in varpr28j. Most computation time ends up being spent 
in the evaluation of the right side of Eq. (|15|l because of 
the triple nested loop required. This is why the access 
time of H(; ; ) and K(\; ) is so important. 

It is preferable to carry out the actual numerical inte- 
gration of Eq. H15fl in terms of the variables wi m (r) , such 
that at m (r) = r wi m {r). This is because these variables 
are constant both at large and small r. 

It is not possible to handle the situation where the 
permittivity inside the NP is a purely negative number. 
That is because smoothing at the boundary would result 
in a point where the permittivity is zero, so that the 
logarithm is singular. However, if the permittivity inside 
is complex valued then this problem is avoided. 



X. NUMERICAL METHODS 

In practice, the expansion for V in Eq. has a limited 
value of £, here denoted as £ m axa- In a similar way, the 
expansion for log(e r ) in Eq. JSJ is limited to a finite value 
of £, denoted l maxc - 

It was found that £ ma xc must be chosen jointly with 
imaxa- Best results are obtained when £ ma xc is less than 
imaxa- Further increasing £ maxa while holding £ ma xc con- 
stant does not improve the results. 

It is desirable to make careful comparisons with the 
exact solution for a homogeneous dielectric ellipsoid in a 
homogeneous matrix. However, in such a situation the 
permittivity is a singular function which is not suitable 
for integration. To get around this problem, the dielec- 
tric ellipsoid is approximated by a smoothed permittivity. 
The smoothing must be carried out over a short distance 
(surface thickness) in order to accurately approximate 



XI. CHECKS OF CORRECTNESS 

The three methods for determining the dipole moment 
given as Eqs. (|4*5|) . (J53}, and (JS7J, do not give identical 
results when finite £ cutoffs ( £ m axa and £ m axc) are used. 
However, mutual convergence of the three methods to 
calculate dipole moment is an indication of convergence 
as the two I cutoffs are increased. 

Given the length of this paper, we have chosen not to 
present specific numerical results. But we can mention 
the numerical tests that we have done as a check of cor- 
rectness of the equations presented here. We compared 
our numerical results to exact ones for the (1) dielectric 
sphere, with positive and complex valued permittivity 
(2) prolate spheroid (3) oblate spheroid and (4) dielec- 
tric sphere with its center not located at the origin of 
coordinates. 

All of these cases have an exact solution that comes 
from Eq. (8.10) on page 41 of Landau's book. 1 - This is 
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for the dipole moment of a homogeneous ellipsoid (with 
semiaxes a, b, and c) in a vacuum with an asymptotically 
uniform electric field along the z-axis. When translated 
from cgs to SI units and into our notation it becomes: 



p z = (volume) Ei- } 



1 + (e* 



l)nW 



(58) 



where is the z-axis depolarization factor which is 
equal to 1/3 for a sphere and "volume" is (47r/3)a6c. 
ti nr is the relative permittivity of the ellipsoid. Using 
our Eq. (|A7|) (in the Appendix) the formula for the case 
of a homogeneous dielectric ellipsoid embedded in a ho- 
mogeneous dielectric is: 



p z — (volume) Ey, 



e (e v 







+ (Cj 



(59) 



For an oblate spheroid where a = b and c < a, the 
eccentricity is e = y a? jb 2 — 1. 



>(*) 



1 



-(e — arctan(e)) 



= n ^ = -(1 -n (z) ) 



(60) 



(61) 



There is an increase in numerical error as the condition 
for dipole surface plasmon resonance is approached. This 
is because -Fioio approaches zero after many positive and 
negative contributions so that numerical errors become 
predominant when I cutoffs are used. 

In order to see the pattern of convergence as £ m axc and 
Pmaxa increase, some series of values have been calculated 
to see if the percentage error goes to zero. The percentage 
error does decrease monotonically. 
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APPENDIX A: HOMOGENEOUS OBJECTS 

We present a fourth method to calculate the dipole 
moment, however this one is specialized to the situation 
of a homogeneous dielectric embedded in a homogeneous 
dielectric. There is no bound charge density pb at interior 
points of a homogeneous dielectric, but bound surface 
charge Ob can exist at dielectric surfaces. To see why, 
note that P = (e r — l))/e r )D, and V • D = pf where pf 
is the free charge density. In a region where e r is constant 
and there is no free charge V • D = so that V • P = 
as well. Thus pb = 0. 



In this case, the volume integral for dipole moment 
p = J rpbd 3 r can be replaced with the surface integral 



P = J Yo- b dA (Al) 
The bound surface charge density is given by 



Pn 



(A2) 



where P ou tn is the normal (outward pointing) component 
of the polarization just outside the surface, while Pi nn is 
the normal (outward pointing) component of the polar- 
ization just inside the surface. This can be related to the 
electric field at the surface since P = e Q (e r — 1)E. Thus: 



0~b e [(€ r i n l)£^ nn (trout ^-)E 



(A3) 



where Ei nn is the normal component of the electric field 
just inside the surface. This is equivalent to: 

0~b Ei nn D ou f n -\- €o(Eoutn Ei nn ) (A4) 

The absence of free charge means that Di nn = D outn . 
Therefore 



o~b e D (E OU f n Ei nn ) 



(A5) 



This provides a formula for dipole moment in terms of 
the electric field: 



p = e a r(E 0Ut - E m ) • dA 



(A6) 



The electric field at great distances is Ei„ c . It can also 
be noted that the electric field in and near the object in 
this situation depends only on ei n /e ou t. To see why, note 
that V • D = 0. Therefore, V • (e(r)E(r)) = 0. Con- 
sequently, a given electric field still solves the boundary 
value problem if both e m and e ou t are multiplied by the 
same constant. 

Because of Eq. (|A6ll , the same also holds true for the 
dipole moment: 



P — /(Ej nc , — — ) 

eout 



(A7) 



where /() is some (non-scalar) function that also depends 
on the size, shape, and orientation of the object. 

Thus, if the dipole moment is theoretically calculable 
for the situation of a homogeneous object sitting in a 
vacuum, Eq. (| AT|) can be used to find the dipole moment 
when permittivities inside and outside are multiplied by 
the same constant. The dipole moment will be the same, 
in fact. 
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